Quantum depletion of collapsing Bose-Einstein condensates 
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We perform the first numerical three-dimensional studies of quantum field effects in the Bosenova 
experiment on collapsing condensates by E. Donley et al. [Nature 415, 39 (2002)] using the exact 
experimental geometry. In a stochastic truncated Wigner simulation of the collapse, the collapse 
times are larger than the experimentally measured values. We find that a finite temperature initial 
state leads to an increased creation rate of uncondensed atoms, but not to a reduction of the collapse 
time. A comparison of the time-dependent Hartree-Fock-Bogoliubov and Wigner methods for the 
more tractable spherical trap shows excellent agreement between the uncondensed populations. We 
conclude that the discrepancy between the experimental and theoretical values of the collapse time 
cannot be explained by Gaussian quantum fluctuations or finite temperature effects. 
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I. INTRODUCTION 

Experimental progress in dilute gas Bose-Einstein con- 
densates (BECs) has recently allowed increasingly de- 
tailed studies of the quantum nature of the atomic 
field [3, [E IE El- This has been accompanied by ad- 
vances in the numerical treatment of many-body quan- 
tum field theory applied to BEC dynamics, most notably 
in a better understanding of phase space methods [5| and 
Hartree-Fock-Bogoliubov theory Q. Therefore, experi- 
ments in which the physics is sufficiently straight-forward 
that quantitative agreement with many-body quantum 
theory can be expected are especially appealing. 

In this paper we extend our previous analyses 0, [1] of 
one of these: the JILA Bosenova experiment of E. Don- 
ley et al. [9(, in which 85 Rb BECs were made to col- 
lapse by switching their atomic interactions to attrac- 
tive. Many interesting phenomena associated with the 
collapse, like bursts and jets, have attracted widespread 
attention. These have been understood qualitatively us- 
ing a variety of models^, & EE El EE EE El El EE 
Il7l EE EE EE HE HE • However precise quantitative 
agreement, of the kind sought in this paper, has not been 
fully achieved. 

Here we are concerned with a quantitative description 
of the most basic aspect of the collapse experiment: the 
time to initiation of the collapse, ^collapse- The abrupt on- 
set of atom losses in the experiment allows a precise mea- 
surement of this time. It has previously been shown that 
the Gross-Pitaevskii (GP) theory substantially overesti- 
mates these collapse times [7] . However quantum correc- 
tions in the framework of time-dependent Hartree-Fock- 
Bogoliubov (HFB) theory were shown not to accelerate 
the collapse of a BEC in a spherically symmetric trap Q . 



Here we investigate the collapse in a cigar shaped trap, 
exactly as in the experiment. Our simulations use the 
stochastic truncated Wigner approximation (TWA), with 
the experimental parameters. We find that the inclusion 
of quantum effects docs not yield results in agreement 
with the experiment. Therefore, for the Bosenova exper- 
iment, despite an excellent qualitative understanding, we 
do not yet have quantitatively precise theoretical models 
for even the simplest aspects. The implications of this 
are discussed in the conclusion. 

We also compare the TWA with the HFB formalism for 
a spherically symmetric trap. We find that the quantum 
depletion predicted by both methods agrees very well. 
As both methods independently confirm an excited state 
population insufficient to accelerate the collapse, we can 
rule out zero temperature depletion as a mechanism for 
collapse acceleration. 

Moreover, using both quantum field methods we show 
that the initial presence of a thermal cloud increases the 
production rate of uncondensed atoms. This results in a 
reduction of the condensate population just before col- 
lapse, which in principle could appear as a slightly ac- 
celerated collapse. However, as we explain in section fVTl 
this effect would be difficult to detect in an experiment. 
Irrespective of this, even for temperatures about three 
times higher than the experimentally measured temper- 
ature of T = 3 nK, the acceleration of the collapse due 
to depletion is insufficient to bring the theoretical and 
experimental results into agreement. 

In summary, we present a careful quantitative study of 
the best characterized experimental data in the Bosen- 
ova experiment: the collapse time. We find that not only 
the GP model, but also the HFB and TWA theories fail 
to explain the collapse times. Further experimental and 
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theoretical work should resolve this unsatisfactory situa- 
tion. 

This paper is organized as follows: In section [IT] we 
give a brief overview of the Bosenova experiment. Sec- 
tion [HT] contains the theoretical background of the TWA 
and HFB methods for the quantum field description of 
BECs. This is followed in section Hvl bv a discussion of 
the numerical limitations of the TWA. Sections Wj and 
IVII contain the results of the TWA for the collapse of a 
BEC in a cigar trap. Finally in section WW we report on 
the comparison between HFB and TWA in a spherically 
symmetric case. 



II. THE BOSENOVA EXPERIMENT 

In the experiment [§], a stable 85 Rb condensate was 
prepared with scattering length a s = using a Feshbach 
resonance, before a s was switched to a negative (attrac- 
tive) value a s = a co ii a pso- The resulting collapsing con- 
densate was observed to lose atoms until the atom num- 
ber was reduced to about the critical value below which 
a stable condensate can exist @. Usually the remnant 
atom number was found to be slightly greater than the 
critical value, a puzzle which has only recently been ex- 
plained, with the help of a new experiment [241 ] . by the 
formation of multiple bright solitons. 

The onset of atom number reduction is quite sudden. 
After the change in scattering length a few milliseconds 
of very little loss is observed. This is followed by a rapid 
decay of condensate population (within ~ 0.5 ms) af- 
ter which the condensate stabilizes again. This behav- 
ior results from the scaling of the loss rate with the 
cube of the density, the peak value of which rises as 
l/(^coiiapse — t) near the collapse point [lGj]. The sud- 
den onset of atom loss allows a precise definition of the 
collapse time i C oiiapse, as the time after initiation of the 
collapse (a s — > a C oiiapse) up to which atom loss remains 
negligible. In this paper we focus our attention primarily 
on a case with a co n apsc = — lOcio, where ao is the Bohr ra- 
dius. For this case the experimentally measured t C oiiapse 
is (6± 1) ms |Sj, |25j, while Gross-Pitaevskii studies found 
it to be about 10 ms 0]. 

A quantitative result of the experiment is the depen- 
dence of tcoiiapse on the magnitude of the attractive inter- 
action, parametrized by the (negative) scattering length 
^collapse- These measurements are performed starting 
with ./Vini = 6000 atoms in an ideal gas state, i.e. the 
interaction between atoms is tuned to zero. The t C oiiapse 
data points presented in [9( have undergone one revi- 
sion of their a co iiapsc values by a factor of 1.166(8) due 
to a more precisely determined background scattering 
length [25|. 

Other experimental features like the bursts and jets 
mentioned in the introduction have also been measured in 
great detail, tempting quantitative explanation. However 
the HFB and TWA methods used in this paper become 
impractical soon after the initiation of collapse, and are 



therefore unsuitable for analysis of the full collapse,even 
if they correctly modelled the collapse times. Refs. 
review the state of theoretical studies of the Bosenova 
experiment. 



III. QUANTUM FIELD MODELS OF A 
HARMONICALLY TRAPPED BEC 

The Hamiltonian for a Bose gas of interacting atoms 
in an external trapping potential V(x) is given by: 

H = J d 3 x * t (x)//o*(x) 

+ ^ j d 3 x *t( x )^rt( x )$( x )4r( x ) j (1) 

£o = -^V 2 + U(x). 

Here ^(x) (^'(x)) is the field operator that annihi- 
lates (creates) a boson at position x, m is the atomic 
mass and Uq — \-Khra s jm is the interaction strength 
with the s-wave scattering length a s . In the following 
we useparameters corresponding to the collapse experi- 
ment [![. The 85 Rb atoms with m = 1.41 x 10~ 25 kg are 
confined in a cigar shaped cylindrically symmetric trap 
V(x) = m(a> 2 _r 2 + w 2 z 2 )/2, where the trap frequencies 
arewi = 17.5 x 2ir Hz and ui z = 6.8x27rHz. For compar- 
ative purposes, we additionally consider a spherical trap 
with the geometric mean frequency Q — {to 2 1 uj z ) 1 /' i = 
12.8 x 2tt Hz in section IVHl 

In the actual BEC collapse atom losses due to three- 
body recombination play a crucial role. These losses are 
taken into account in the master equation for the time 
evolution of the system's density operator p 26] : 

-* t (x) 3 *(x) 3 p-p* t (x) 3 #(x) 3 ). (2) 

The three-body loss constant K3 = 1 x 10 ms is 
chosen as in |7fl. K% is not very well constrained exper- 
imentally, but in simulations its value can be varied by 
factors of 10 to 100 without significantly affecting the col- 
lapse time ■ The reason for this is that the three-body 
loss acts only as a diagnostic for a rapid increase in den- 
sity at the point of collapse. In the time leading to this 
increase, the density of the contracting BEC remains low 
enough for three-body loss to play no role in the dynam- 
ics. Only at and after the actual time of collapse does 
the precise value of K 3 become relevant. 

In the following subsections we describe two different 
approaches to finding approximate solutions of the quan- 
tum evolution given by Eq. ([2"|). 
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A. Truncated Wigner method 

To obtain the time evolution of p we may represent p 
in a suitable phase-space ||. In this case we make use 
of the Wigner representation. We define the multimode 
Wigner function: 



W({a k }, {a* k }) = II (7 \d 2 (3 h \ cxp ( ]T /?, 



Tr[exp(^/3,aj-/J 



3 a 3 ) P 



(3) 



Here aj (%) creates (annihilates) atoms in the jth single 
particle mode. These may be eigenstates of the harmonic 
oscillator or position eigenstates on a discrete grid. Using 
this representation it is possible to obtain the evolution 
of ^({afe}, {a£}) from Eq. If the resulting equation 
is truncated by neglecting derivatives higher than second 
order with respect to the a^, it takes the form of a Fokker- 
Planck equation (FPE), which can then be mapped onto 
a stochastic differential equation (SDE) Q . The validity 
criteria for the truncation and details of the procedure 
are discussed in Refs. [13, [H, HsJ. 

Other theoretical approaches, such as the positive-P 
or gauge-P phase space methods [30|, HH , can include the 
full quantum evolution of the s-wave scattering physics, 
but still necessitate an approximate description of three- 
body losses. We have implemented both these methods 
for a one dimensional, and also for a spherically symmet- 
ric, collapse scenario and found them to be numerically 
intractable. 

Compared to the situation without loss [13] the inclu- 
sion of three-body recombination in the master equation 
results in additional terms in the stochastic differential 
equation. These have been thoroughly treated in [32| . 
Drawing on these pre vious studies, we can write down 
the simplest SDE [33[ which describes a trapped BEC 
with three-body loss: 



#(x) 



2m 



V 2 + V(x) + C/ o |0(x)| 2 Wx)di 



^|0(x)|V(x)dt 



^Wx)| a deOr,t). (4) 



Details regarding the construction of the dynamical noise 
term d£(x, t) are given in appendix [X] By setting 
d£(x,t) = in Eq. ^ we can recover the usual Gross- 
Pitaevskii equation including three-body loss [7|. 

The initial state of the stochastic wavefunction <f>(x) 
has to be chosen such that it represents the Wigner func- 
tion of an initial coherent state BEC. At zero temper- 
ature this is achieved by the addition of initial vacuum 
noise ry(x): 



(x,t = 0) = Vo(x) 



7r (x) - 



(5) 



Here V'oO*) denotes the oscillator ground state, which 
is appropriate since the starting point of the experi- 
ment is a non-interacting BEC. r/(x) is a Gaussian dis- 
tributed complex random function that fulfills the condi- 
tions 7j(x)t/(x') = and r)(x)*r)(x') — <5(x — x'), where / 
denotes the stochastic average of /. 

The truncation of higher order derivatives in the FPE 
is safely applicable when all modes in the problem are 
highly occupied [28[. If the three-dimensional collapse 
scenario is numerically solved on a spatial grid as in [3] , 
6000 atoms are spread out over ~4x 10 6 position space 
modes (i.e. grid points). Thus in the position basis the 
mode occupation criterion cannot be fulfilled. Also the 
addition of the initial noise as in Eq. ((SJ becomes prob- 
lematic, leading to aliasing effects at the edge of the com- 
putational grid. These can be overcome in periodic situ- 
ations as described in [34| . but would be persistent in a 
harmonic trap. 

Here we choose a powerful method to solve the Gross- 
Pitaevskii equation in the oscillator (energy) basis in- 
stead. The method was presented in 13511 and applied 
to BECs at finite temperature in [3(| |37| |. In our work 
we extend the formalism in order to solve Eq. (H|). The 
stochastic field </>(x) in Eq. ([4} can be expanded in terms 
of eigenstates if of the ID harmonic oscillator: 

C-lm.n (t)<pi(x)<p m (y)ip n (z), (6) 

{l,m,n}£C 



which fulfill: 



h 2 d 2 1 



Similarly if m {y) an d ^ n {z) solve the oscillator equation 
in the y- and z-dimensions. The summation in Eq. ([6]) is 
restricted to all modes with total energy below a certain 
cutoff E cut : 



C = {l,m,n: e; + e m + e n < E cut } . 



(8) 



The values for i? cu t are given in units of hiu± (to) for 
the cylindrical (spherical) case. Further details of the 
approach are given in appendix [XJ 

The stochastic equation dU) has to be solved for many 
realizations (trajectories) [38|]. If the Wigner represen- 
tation is used, symmetrized quantum averages can then 
be determined from averages over all trajectories Q. In 
what follows (/) denotes the quantum expectation value 
of /. We obtain the condensed (coherent) and uncon- 
densed (incoherent) populations for each oscillator mode 
from: 



Aran 
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(9) 



(10) 
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The i>imn (*f„) are field operators that annihilate 
(create) an atom in the mode with quantum numbers 
l,m,n. The numbers of condensed and uncondensed 
atoms (-/V con d, iV unc ) are then obtained by summing the 
populations in all the modes. The total atom number is 
given by iV tot = iV con d + N unc . 

A more rigorous definition of the condensate com- 
ponent of the stochastic field is given by the Penrose- 
Onsager criterion [39[. Exemplary applications of this 
method can be found in [36|, [40] . To employ the criterion 
we would need to average and subsequently diagonalizc 
the one-body density matrix of size N mo dea x ^modesi 
which however is not feasible in our case as will be ex- 
plained in section [Tvl 

Finally we point out that the mode occupation cri- 
terion of [28| can be slightly relaxed to the require- 
ment that the noise density defined by <5 c (x, x) = 
£{J,m,n}eC \ ( Pi( x ) l Prn(y)'Pn(z)\ is smaller than the con- 
densate density n c within the volume where the latter 
is significant [2^, |4l| . This criterion is basis dependent 
and for the Bosenova problem it can be fulfilled in the 
oscillator basis but not in the position basis. 



B. Time-dependent Hartree-Fock-Bogoliubov 
approach 

A different method to go beyond the mean-field the- 
ory is to derive the Heisenberg equation for the field 
operator ^ a (x, i), and subsequently decompose ^(x, t) 
into a condensate part <p a (x, t) and quantum fluctua- 
tions x(x, t), such that *P a — <Pa + X and (^f a ) = <f> a . 
The quantum fluctuations can be described in terms of 
their lowest order correlation functions: the normal den- 
sity Gjv(x, x') = (x^( x ')x( x )) an d anomalous density 
Ga(x,x') = (x(x')x(x)). The derivation of the dynam- 
ical equation for the condensate contains a factorization 
of the expectation values in accordance with Wick's the- 
orem g3, e.g. (&JfJf a ) =2<#J *„><*„} + <#t}<# * a >. 

This implies the assumption that the system is in a 
Gaussian quantum state (i.e. a coherent state or even 
a squeezed state). We obtain as dynamical equation for 
the condensate: 



To obtain the time evolution of the condensate we have to 
supplement Eq. (fTTj) by evolution equations for Gn and 
Ga, which are listed in appendix iBl We refer to Ref. [§| 
for further details regarding a method to solve the set of 
coupled equations in a harmonic trap. These equations 
require a renormalisation of the coupling strength Uo due 
to the momentum cutoff K — it /Ax of the numerical grid 
used in the simulations [43j]. One must distinguish the 
physical interaction strength U and the parameter Uq in 
the Hamiltonian, which are related by: 



Uq = 



u 



1-aU 1 



mK 



(13) 



A similar renormalisation issue arises in the truncated 
Wigner method where the same prescription has to be 
used to relate the numerical coupling to the physical in- 
teraction strength [28| . In the interaction strength regime 
of interest for this paper the difference between U and Uq 
is negligible. A careful renormalisation is hence unneces- 
sary and we can directly employ Uq = U = Airh 2 a s /m in 
our simulations. 

The cutoff is determined by the spatial lattice spacing 
in the HFB case, and by the energy cutoff E cut in the 
TWA case, after equating the free particle kinetic energy 
to the energy cutoff: E cut = h 2 K 2 / (2m). The cutoffs are 
chosen to ensure numerical accuracy of the simulations, 
and in particular that the results of interest are invariant 
with respect to changes in the cutoffs. For our HFB 
simulations the highest cutoff is K = 1.3 x 10 7 m -1 , and 
\aU\ — 4.5 x 10 -3 , corresponding to less than one-half 
percent renormalization. For our TWA simulations the 
maximum cutoff is E cut « 50 and therefore K = 3.8 x 10 6 
m _1 , and \aU\ = 1.3 x 10~ 3 , again corresponding to 
negligible renormalization. 

We note that the three-body loss term in the HFB 
formalism, Eq. (|12[) . only incorporates loss processes 
among condensate particles, whereas the implementation 
in Eq. ((4]) of the TWA contains loss also for the uncon- 
densed fraction. The comparison between the methods 
in this paper is done for very small uncondensed popula- 
tion and small total losses, so that this difference can be 
neglected. 



■ fc <90 a (x) 



1 2m 



+ 2UoM*)Gn (x, x) + (/oC(x)G A (x, x) . (11) 

Here Gat (x, x) represents the density of uncondensed 
atoms. Hence this modified Gross-Pitaevskii equation 
contains the interaction of the uncondensed component 
with the mean-field. 

We also phenomenologically model three-body loss 
from the condensate, by adding the following term to 
equation (fTT|) : 



qK 3 \u*)\ 4 Mx). 



(12) 



IV. NUMERICAL CONSTRAINTS 

The aim of this section is to clarify the basis of the 
conclusions drawn from our simulations. 

In this article we present solutions of Eq. ^ , model- 
ing the Bosenova experiment without any significant free 
parameters, using three different levels of approximation 
of the truncated Wigner method: 

• GP evolution: Gross-Pitaevskii evolution only. In 
this case both the noise on the initial state and 
dynamical noise are omitted (rj = 0, g?£ = 0). 

• TWA with initial noise: Truncated Wigner evolu- 
tion without dynamical noise (tf£ = 0). 
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• TWA with dynamical noise: Complete truncated 
Wigner evolution (77 ^ 0, d£ 7^ 0). 

The reasons for studying the GP evolution are two- 
fold. Firstly it aids the determination of the required 
number of oscillator modes by comparison with the es- 
tablished position space results [7(. Secondly it allows 
us to quantify the differences between the classical and 
quantum field results. 
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FIG. 1: GP evolution only. Atom number N con £ in the con- 
densate during collapse with a co iiapsc = — 10ao- (•) Results 
obtained on a spatial grid 0. Thick lines correspond to so- 
lutions of Eq. Q in the energy basis. Inset: A table with 
the number of modes for different E cu t and a legend. As 
S cu t increases, the atom number curves approach the correct 
position basis solution. 

To determine the required mode numbers, the GP 
equation is solved in the harmonic oscillator basis in or- 
der to reproduce the atom number curve of Ref. 0- In 
doing so, we encountered a limitation of the oscillator- 
basis: due to the extremely narrow peak of the conden- 
sate wavefunction at the collapse time [l4j], numerically 
accurate simulations beyond this point require a very 
large number of modes, > 10 6 . Simulating the conden- 
sate evolution much beyond the collapse time is therefore 
not feasible [44[ ■ Fig. Q] shows the number of atoms re- 
maining in the condensate for different numbers of oscil- 
lator modes employed. In the case of 4 x 10 5 modes, the 
result appears close to convergence against the solution 
of the GP obtained on a spatial grid [3] • However, we can 
conclude from the evolution of the peak densities that a 
cutoff of at least E cut > 150, corresponding to about 
1.5 x I0 6 modes, would be required to evolve through the 
collapse. Nevertheless, we find that the evolution until 
~ 8 ms can be accurately represented with a basis-size 
that is computationally tractable in the stochastic multi- 
trajectory treatment (~ 5 x 10 4 modes, E cut = 50). With 
this mode-number, the validity criterion of [2!| is safely 
fulfilled. 



interaction. Since the interaction between the condensed 
and uncondensed components is twice as strong as the 
self-interaction of the condensate (compare Eq. (jTTJ) ) , a 
sufficiently strong quantum depletion could possibly yield 
a quicker collapse than a pure BEC. 

Here we present results for a collapse with a co n ap sc = 
— IOao, for which the experimentally measured i C oiiapse is 
(6 ± 1) ms 0, [25| . We find that only very few uncon- 
densed atoms are created prior to the actual collapse, see 
Fig. [5] (a) . This result qualitatively agrees with our pre- 
vious studies of a spherically symmetric geometry with 
the HFB method Due to the very small depletion in 
the initial stage, our quantum treatment does not show 
an acceleration of the collapse compared to the GP evo- 
lution. We deduce this from the identical peak densities 
in Fig. [5] (b). As discussed in the previous section, de- 
spite the numerical limitations we can thus conclude that 
the inclusion of zero temperature quantum depletion does 
not result in agreement between theory and experiment. 

Fig. [2] (a) shows that the dynamical noise contributes 
significantly to the evolution of the uncondensed atom 
number. Fig. 0] (b) confirms that dynamical noise is nec- 
essary for exact agreement with the HFB approach. 




FIG. 2: Initial 8 ms of evolution after the change in scat- 
tering length from to a co n a p SO = -10a . (a) iVunc for 
the solutions with dynamical noise (solid) and initial noise 
only (dashed), (b) Condensate peak density for GP evolution 
(solid) and TWA with dynamical noise (dashed). The result 
is unchanged, hence no acceleration of the collapse occurred 
before 8 ms, which exceeds the experimental collapse time of 
(6 ± 1) ms. In both (a) and (b) E cut = 50. 

We have checked that these results are qualitatively 
unchanged for other scenarios, i.e. where the attractive 
interaction is changed to a co ii a psc = — 6ao and a co iiapsc = 
— 25ao- In the latter case, we have observed a larger 
production of uncondensed atoms just before collapse, 
which also did not significantly accelerate the collapse. 



VI. INCLUSION OF A THERMAL CLOUD 



V. COLLAPSE OF A CIGAR SHAPED BEC 

If Eq. ((4]) is solved in the truncated Wigner formal- 
ism, the condensed and the uncondensed fractions of the 
atomic gas can be distinguished. During the collapse 
population is transferred between the fractions due to the 



As a next step towards even more accurate modelling 
of the experiment we include initial thermal population 
in the uncondensed modes. We will show that, if tem- 
perature effects are taken into account, the precise results 
for the measured collapse times might depend on whether 
N con d or JVtot is measured. 
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In the oscillator basis representation, the initial state 
for nonzero temperature [28| is 



Clr. 



■00, 



ran I Vlrnn 



2 tanh 



2k B T 



-1/2 



(14) 



Here e; m „ = e; + e TO + e« are the energies of the os- 
cillator modes with quantum numbers I, m, n. T is the 
temperature of the cloud and i/'o.Zmn are the correspond- 
ing expansion coefficients of the GP ground state. The 
Vimn are complex Gaussian noises obeying r]f r)ii m ' n ' = 
Sii'S mm '8 nn i. Although the temperature in the Bosenova 
experiment was measured to be 3 nK, in Fig. [3] we present 
our results for a few different temperatures: T = 0, 
T = 3 nK, T = 5.3 nK and T = 8 nK. We plot the final 
2 ms of simulated time. This corresponds to the collapse 
stage and exceeds the time of ~ 8 ms for which we can 
employ sufficiently many modes to ensure a reliable sim- 
ulation. Nonetheless we would like to draw some qualita- 
tive conclusions from these simulations of "BEC collapse 
in a restricted mode space" . Firstly, if the collapse time 
was deduced from the condensate atom number alone, 
which is shown in Fig. [3] (a), it appears to be shorter for 
increased temperature. The reduction by ~ 0.75 ms for 
the experimental temperature of 3 nK is however by far 
not enough to reach agreement with the experimental col- 
lapse time of (6 ± 1) ms. Secondly, the reduction in con- 
densate atom number just before collapse, compared to 
the GP dynamics, results from stimulated transitions to 
uncondensed modes rather than an increased total atom 
loss, which can be deduced from an inspection of N unc 
and -/Vtot- These qualitative features are independent of 
the value of E cut used in the simulations. However we 
point out that the quantitative details of the evolution of 
condensed and uncondensed fractions for the times pre- 
sented in Fig. [3] depend on E cut . Fig. [3] (b) shows that 
the acceleration of collapse is much smaller, if only the 
total atom number is taken into account. 
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FIG. 3: Slight acceleration of the collapse due to initial ther- 
mal atoms, (a) Time evolution of AT con( j around the collapse, 
(b) Change in total number AN = iV to t (t) - N to t (0) for the 
same period of the evolution. The sampling errors of all these 
results are less than 1% [ill ]. £" C ut = 50. 

In the experiment [9( , the measured atom number was 
deduced from counting atoms in the central region of 
the trap. The increased population of the uncondensed 



modes just before collapse occupies the same spatial re- 
gion as the condensate. Since the experimental method 
did not distinguish between condensed and uncondensed 
atoms, we conclude that the experiment probably did not 
capture the above described temperature effect. We note 
that, for T / 0, the evolution of the total number of 
atoms is only marginally changed compared to the GP 
results. 



VII. HFB VS. WIGNER: COLLAPSE OF A 
SPHERICAL BEC 

In this section we compare the two different quantum 
field models of BECs used in this paper. Both rely on ap- 
proximations to achieve a numerically tractable descrip- 
tion of the quantum evolution. The formalism of each 
method and the approximations involved differ greatly, 
as outlined in section IIIII The quantitative agreement 
between the evolution of the uncondensed fraction in 
both methods that we present in this section gives thus 
a strong indication of the validity of our results. 

As has been described in Ref . [8| , Hartree-Fock Bogoli- 
ubov simulations are not feasible in the case of the cylin- 
drically symmetric experimental situation, as the correla- 
tion functions G a and Gn become then five dimensional. 
Therefore in our earlier work [8j, we have used the HFB 
to investigate a collapse in a spherically symmetric trap. 
For the same reasons, we compare here the TWA and 
HFB methods for a range of temperatures of the initial 
thermal cloud in a spherical geometry. Finite tempera- 
ture is taken into account in the HFB approach by using 
correlation functions corresponding to a thermal popula- 
tion of oscillator states initially: 



G N (x,x') = ]T 

It 

G A (x,x')=0 



1 



, expf^V^) - 1 

Imn ^ y k B l ' 



^*mn( x ')¥W( x )> 



(15) 
(16) 



Here ipj mn (x) = tpi(x)ip m (y)tp n (z). For the finite tem- 
perature TWA simulations we have used E cu t = 40. 

We find excellent agreement between the uncondensed 
atom numbers predicted by HFB and TWA for the initial 
5 ms of the collapse [46[ and a range of temperatures, 
as shown in Fig. |4j The higher the temperature, the 
larger the initial uncondensed population iV unc (0), which 
causes more stimulated transitions to the uncondensed 
fraction. As reported in |47| an increase in temperature 
of the initial state also requires more trajectories for the 
sampling error to be satisfactorily small. 

The TWA and HFB in the spherical case both qual- 
itatively confirm the results presented in section I VII for 
cylindrical geometry: higher temperature thermal clouds 
yield an increased creation of uncondensed particles just 
before collapse. This appears like an accelerated collapse 
on the curve for iV con( j. In contrast, inspection of the 
total atom number shows almost no acceleration. 
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Validity timescale: For the period we consider here, the 
approximations involved in both methods are justified 
and therefore a comparison is useful. It is known that 
the truncation in the TWA and the factorization of cor- 
relation functions in the HFB are valid for short times 
only, but both methods suffer from a lack of quantita- 
tive knowledge about this timescale in the general case. 
For the situation of a BEC in an optical lattice within 
the tight binding (Bose-Hubbard) regime, the validity 
timescales for TWA and HFB have been shown to coin- 
cide. They are given by t <C J/U, where J and U are the 
Bose-Hubbard hopping strength and on-site interaction 
respectively (48l. l49jj. 




2 3 
t(ms) 



FIG. 4: Increase in uncondensed atom number N U nc{t) — 
A^unc^O) during the first stages of a spherical collapse with 
a C ollapsc = -12oo for the TWA and HFB. (a) The TWA re- 
sults with dynamical noise (dashed) are in excellent agreement 
with the HFB result (solid) for T = 0,3, 5.3 nK. Dotted lines 
indicate the sampling error. Numerical parameters are given 
in the footnote [50]. (b) Close-up of the result for T = 0. 
As mentioned in section [V] the TWA with dynamical noise 
(dashed) agrees better with the HFB (solid) than the result 
of the TWA with initial noise terms only (dot-dashed). 



VIII. CONCLUSIONS 

Our three dimensional simulations of the Boscnova ex- 
periment on collapsing BECs 9] have shown a moderate 
acceleration of the collapse if an initial thermal cloud 
is taken into account, although the effect is not large 
enough to quantitatively account for the discrepancy be- 
tween experimental and theoretical collapse times. The 
predictions of Hartree-Fock Bogoliubov and truncated 
Wigner theories for collapse in a spherically symmetric 
case agree very well with each other. However in a cigar 
shaped trap, where only the truncated Wigner method is 
feasible, the theory disagrees with the experiment. The 
origin of this discrepancy is unknown. 

Close to a Feshbach resonance molecular effects can 
become important. However during the sequence of the 
Bosenova experiment, the magnetic field stays clear of 
the resonance value [9(. Hence molecules play a minor 
role in this experiment, as already argued before [1, 
Also a possible breakdown of the s-wave approximation 
and strong effects due to inelastic collisions do not oc- 
cur prior to the actual collapse and hence cannot affect 
the collapse time. As other options are ruled out, we 
conjecture that the collapse time discrepancy arises from 
quantum correlations not captured by our descriptions. 
Although a Gaussian initial quantum state is physically 
reasonable, high order correlations can be created by the 
interactions and could become significant during collapse. 
These are not captured by the methods employed here. 

To understand the discrepancy better, it would be de- 
sirable to conduct an experiment with the aim of mea- 
suring the collapse times and correlation functions for a 
larger range of scenarios and with even higher precision. 



APPENDIX A: PROJECTED 
GROSS-PITAEVSKII EQUATION IN THE 
ENERGY BASIS 



Numerical performance: We find that both methods for 
the quantum field treatment of Bose gases, HFB and 
TWA, agree in a direct comparison of the uncondensed 
atom number evolution during a BEC collapse. The 
TWA is advantageous for spatially asymmetric problems, 
as the increasing dimensionality of the correlation func- 
tions renders the HFB method numerically intractable. 
However, in a spherically symmetric case the HFB is ad- 
vantageous. This is because the correlation functions 
Gjv and G a are only three dimensional due to the spa- 
tial symmetry. Meanwhile the truncated Wigner evolu- 
tion has to be always modelled in full three dimensions 
as the quantum fluctuations cannot be assumed to be 
spherically symmetric. The quantum evolution in the 
HFB approach is obtained with a single solution of equa- 
tions (fTTjl . (|B2[) and (|B1[) . compared to averaging over 
many realizations of Eq. (j4|) in the TWA. As a result our 
simulations showed that for the spherical case the HFB 
requires vastly shorter CPU times [HTj] . 



Using the expansion ^ the stochastic equation (|3|) 
becomes: 

dcimn = — — (e; + e m + ^o^lmn) dt 



TrGimndt + J ——dH lmn . (Al) 



The F, G and dH are overlap integrals defined by: 





f d 3 x 




^(y)^(z)|0(x)| 2 0(x), 


(A2) 




f d 3 x 




^(y)^(z)|0(x)|V(x), 


(A3) 




1 d 3 x 


<P*( X 


<(y)<(z)|0(x)| 2 de(x), 


(A4) 


de(x) = 


E 


d£i mn (t)ipi(x)ip m (y)(p n (z). 


(A5) 



{l,m,n}eC 



d£imn(t) in Eq. (| A5|) are complex Gaussian noises that 
fulfill d£i mn (t)d& m , n ,(t') = and d£f mn (t)d&< m '„/(i) = 

Sll'S mm 'Snn'dt. 

It has been outlined in Refs. [H, Hi| how these in- 
tegrals can separately be exactly computed on an ap- 
propriately chosen non-equidistant spatial grid. Differ- 
ent spatial grids would however be necessary for the ex- 
act solution of integrals involving different powers of the 
wavefunction. To remain computationally efficient we 
have chosen a grid which allows the exact calculation of 
Eqns. (|A2|) and (|A4[) . We have checked that our results 
are invariant under a variation of the grid used for eval- 
uation of the integrals. 
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ih d -^^l = {[x{x')x{x),H]) = 

))G A (X,X') + 2[/ ( |0a(x)| 2 + |0a(x')| 2 

+ G N (x, x) + G N (x', x') ) G A (x, x') 

+ Uo(Mx) 2 G* N (x, x') + M*') 2 Gn (x, x') 

+ G A (x, x) G* n (x, x') + G A (x, x') Gat (x, x') ) 

+ C/ o (0 a (x) 2 + G A (x,x) )<5^(x - x'). (B2) 



APPENDIX B: HFB EQUATIONS FOR 
CORRELATION FUNCTIONS 



In the time-dependent Hartree-Fock-Bogoliubov ap- 
proach we have to supplement Eq. (|II[) by evolution equa- 
tions for Gn and Ga- These are obtained by deriving the 
Heisenberg equations for the operators x\ x ')xi x ) an d 
x( x ')x( x ) respectively, and factorizing operator products 
as described in Q. This procedure yields: 

ih dG N (g*) =m j mx)tA]) = 

(Hoa(x) - Hoa(x'))G N (x,x') + 2U ( |0 Q (x)| 2 - |0 a (x')| 2 

+ G N (x, x) - Gat (x', x') )G n (x, x') 

+ U (G A (x,x) G* a (x,x) - G* A (x',x) G A (x,x) ) 

+ U {M*) 2 G A (x,x') - ^) 2 G A (x,x')), (Bf) 
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